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ABSTRACT 

One of the most important unresolved issues in gamma-ray burst physics is the origin of the prompt 
gamma-ray spectrum. Its general non-thermal character and the softness in the X-ray band remain 
unexplained. We tackle these issues by performing Monte Carlo simulations of radiation-matter in¬ 
teractions in a scattering dominated photon-lepton plasma. The plasma - initially in equilibrium - is 
driven to non-equilibrium conditions by a sudden energy injection in the lepton population, mimicking 
the effect of a shock wave or the dissipation of magnetic energy. Equilibrium restoration occurs due 
to energy exchange between the photons and leptons. While the initial and final equilibrium spectra 
are thermal, the transitional photon spectra are characterized by non-thermal features such as power- 
law tails, high energy bumps, and multiple components. Such non-thermal features are observed at 
infinity if the dissipation occurs at small to moderate optical depths, and the spectrum is released 
before thermalization is complete. We model the synthetic spectra with a Band function and show 
that the resulting spectral parameters are similar to observations for a frequency range of 2-3 orders of 
magnitude around the peak. In addition, our model predicts correlations between the low-frequency 
photon index and the peak frequency as well as between the low- and high-frequency indices. We 
explore baryon and pair dominated fireballs and reach the conclusion that baryonic hreballs are a 
better model for explaining the observed features of gamma-ray burst spectra. 

Subject headings: gamma-ray burst: general — radiation mechanisms: non-thermal 


1. INTRODUCTION 

The radiation mechanism that produces the bulk of the 
prompt emission of Gamma-Ray Bursts (GRBs) is still 
a matter of open debate (e.g. Mastichiadis & Kazanas 
2009; Medvedev et al. 2009; Ryde & Pe’er 2009; Asano 
et al. 2010; Ghisellini 2010; Lazzati & Begelman 2010; 
Daigne et al. 2011; Massaro & Grindlay 2011; Resmi 
& Zhang 2012; Hascoet et al. 2013; Crumley & Ku¬ 
mar 2013). Among the many proposed possibilities, the 
synchrotron shock model (SSM) and the photospheric 
model (PhM) have recently gathered most of the atten¬ 
tion (Rees & Meszaros 1994; Piran 1999; Lloyd & Pet¬ 
rosian 2000; Meszaros & Rees 2000; Rees & Meszaros 
2005; Giannios 2006; Pe’er et al. 2006; Bosnjak et al. 
2009; Lazzati et al. 2009; Beloborodov 2010; Mizuta 
et al. 2011; Nagakura et al. 2011). Within the SSM, 
the bulk of the prompt radiation is produced by syn¬ 
chrotron from a non-thermal population of electrons gy¬ 
rating around a strong, locally-generated magnetic field. 
The non-thermal leptons are produced either by trans- 
relativistic internal shocks (the SSM proper, Rees & 
Meszaros 1994) or by magnetic reconnection in a Poynt- 
ing flux dominated outflow (e.g. the ICMART model, 
Zhang & Yan 2011). The SSM naturally accounts for the 
broad, non-thermal nature of the spectrum. However, it 
has difficulties in accounting for bursts with particularly 
steep low-frequency slopes (Preece et al. 1998; Ghisellini 
et al. 2000) and has limited predictive power, since the 
radiation properties are tied to poorly constrained quan¬ 
tities such as the lepton’s energy distribution, the ad- 
hoc equipartition parameters, and the ejection history of 
shells from the central engine. 

The PhM does not specify a radiation mechanism, as¬ 


suming instead that the burst radiation is produced in 
the optically thick part of the outflow and advected out, 
its spectrum being the result of the strain between mech¬ 
anisms that tend to bring radiation and plasma in ther¬ 
mal equilibrium and mechanisms that can bring them out 
of balance (e.g., Beloborodov 2013). The PhM has been 
shown to be able to reproduce ensemble properties of the 
CRB population, such as the debated Amati correlation, 
the Golenetskii correlation, and the recently discovered 
correlation between the burst energetics and the Lorentz 
factor of the outflow (Amati et al. 2002; Amati 2006; 
Liang et al. 2010; Fan et al. 2012; Ghirlanda et al. 
2012; Lazzati et al. 2013; Lopez-Camara et al. 2014). 
However, it is not yet understood how the broad-band 
nature of the prompt spectrum, spanning many orders 
of magnitude in frequency, is produced. In a hot, dissi¬ 
pationless flow, only the adiabatic cooling of the plasma 
would work as a mechanism to break equilibrium, and 
the CRB outflow would work as a miniature big bang, 
the entrained radiation maintaining a Planck spectrum. 
In a cold, dissipationless outflow, lepton scattering domi¬ 
nates the radiation-matter interaction producing a Wien 
spectrum (Rybicki & Lightman 1979). Outflows from 
CRB progenitors are, however, far from dissipationless. 
Hydrodynamic outflows are continuously shocked out to 
large radii (Lazzati et al. 2009), and Poynting-dominated 
outflows suffer dissipation through magnetic reconnec¬ 
tion (Giannios & Spruit 2006). Either way, even if ther¬ 
mal equilibrium is reached at some point in the outflow, 
it is likely that such equilibrium is broken by a sudden 
release of energy in the lepton population or altered by 
a slow and continuous (or episodic) injection of energy. 
The effects of such energy injection on the photospheric 
spectrum are profound (e.g., Giannios 2006; Pe’er et al. 
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2006; Beloborodov 2010; Lazzati & Begelman 2010). In 
addition, the interaction between different parts of the 
outflow in a stratified flow alter the thermal spectra into a 
non-thermal, highly polarized spectrum (Ito et al. 2013, 
2014; Lundman et al. 2013). 

In this paper we investigate the evolution of the radi¬ 
ation spectrum following the sudden injection of energy 
in the lepton population of a plasma, assuming that the 
radiation and leptons interact via Compton scattering 
and pair processes. We use a Monte Carlo (MC) method 
that evolves simultaneously the photon and lepton pop¬ 
ulations by performing inelastic scattering between pho¬ 
tons and leptons in both the non-relativistic and the 
relativistic (Klein-Nishina) regimes. The code also ac¬ 
counts for e^e"*" annihilation (pair annihilation hence¬ 
forth) and e^e"*" pair production from photon-photon col¬ 
lisions (pair production henceforth). We focus on tran¬ 
sient features that can be observed if the episode(s) of 
energy injection in the leptons occur at small or moder¬ 
ate optical depths (r < 1000). 

This manuscript is organized as follows. In Section 
2 we describe the physics and the methods of the MC 
code, in Section 3 we show our results and in Section 
4 we discuss the results and compare them to previous 
Hndings. 


2. METHODOLOGY 
2.1. Step 1: Particle Generation 

As a first task, the code generates a user-defined num¬ 
ber of leptons and photons. Their energies follow a distri¬ 
bution that can be either of thermal equilibrium (Wien 
for the photons and Maxwell-Jiittner for the leptons) or 
any other user-specified distribution. After initializing 
the photon and lepton distributions, our code performs 
the following steps iteratively. 

2.2. Step 2: Particle/Process selection 

To initiate either a scattering or a pair event we need 
to select two particle^ - which we obtain by randomly 
selecting a pair from our generated distributions. De¬ 
pending upon the particles selected, Compton scattering 
(if a photon and a lepton is chosen), pair annihilation (if 
an e~ or e"*" is chosen) or pair production (if two photons 
are chosen) is performed or another pair is re-selected if 
any other combination occurs. After the selection, the 
code proceeds with the following calculations: 

1. Incident angle generation (9) using the appropriate 
relativistic scattering rates, under the assumption 
that both leptons and photons are isotropically dis¬ 
tributed. 

2. Lorentz boost to the necessary reference frames 
(details explained in successive sections) from the 
lab frame. 

3. Event probability computations from total cross 
section (a) calculations. 

4. Scattering angle generation from differential cross 
section (^). 

^ Note that here particle can mean both a lepton or a photon. 


5. Lorentz boost from the necessary frame back to the 
lab frame. 


In the following sub-sections we discuss each of the three 
possible processes in detail 


2.2.1. Process 1: Compton Scattering 

As the choice of reference frame is arbitrary, in the lab 
frame we can assume that the lepton is traveling along 
the x-axis and the photon is incident upon the lepton in 
the xy plane without any loss of generality. The angle of 
incidence between the chosen photon-lepton pair is 
generated by a probability distribution P^e- 

^70) OC sin^.Ye(l /^e COS^^e) (1) 

where f3e = Ve/c, is the ratio of lepton speed to the speed 
of light. 

To simulate the scattering event the code Lorentz trans¬ 
forms to the lepton frame (that we call the co-moving 
frame). The probability that the chosen photon-lepton 
pair interacts depends on the incident photon energy in 
the co-moving frame. As Compton scattering becomes 
less efficient at higher energies, photons having energies 
comparable to or greater than the lepton’s rest mass en¬ 
ergy are less likely to scatter. Using Monte Carlo sam¬ 
pling we determine if scattering occurs or not. This is 
done by generating a random number and comparing it 
to the ratio of the Klein-Nishina cross section a^e to the 
Thomson cross section, which we use as a reference value. 
We proceed with the scattering event if > si 

where si is a random number. If the condition is not 
satisfied, the code returns to step 2. If instead the condi¬ 
tion is satisfied and the scattering occurs, the code gen¬ 
erates the polar scattering angle in accordance with 
the Klein-Nishina differential cross-section formula 


da-y 
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2 E'^ 
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( 2 ) 


where ro is the classical radius of an electron, E' and 
Eg are the energies of the incident and scattered pho¬ 
ton respectively (e.g. Bluementhal & Gould 1970, Lon- 
gair 2003 and Rybicki & Lightman 1979). The energy 
transfer equation connecting E' with Eg is the Compton 
equation (e.g. Bluementhal & Gould 1970, Longair 2003 
and Rybicki & Lightman 1979) 
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( 3 ) 


(Note here that 0' is the angle that the scattered pho¬ 
ton makes with the direction of propagation of the inci¬ 
dent photon in the co-moving frame. Hence equations ([2]) 
and ([3|) hold true only in the lepton frame). Finally, the 
azimuthal angle is generated randomly between zero 
and 27r. Thus, we now have the four momenta of the 
scattered particles in the co-moving frame. 


2.2.2. Process 2: Pair Production / Photon Annihilation 

If the particle selection process selects two photons 
then the pair production/photon annihilation channel is 
chosen. The code computes the angle of incidence 
between the chosen photons by using the probability dis¬ 
tribution Pyy\ 

Py-y{9jy) OC SiU (1 — COS9yy). 


( 4 ) 
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To ensure that the photon pair has enough energy to 
lead to a pair production event the code checks the en¬ 
ergy of the photon/s in the zero momentum frame. The 
zero momentum frame photon energy E'^ can be com¬ 
puted given the incident photon energies Ei , E 2 and the 
incident angle as 

E'^ = /^sin(0^^/2). (5) 

(Gould & Schereder 1967). If £1' < meC? the colliding 
photon pair is not energetic energy to produce an e“e+ 
pair, hence the code jumps to step 2 for a new particle 
pair selection. Due to the energy dependence of cross- 
section , even photons exceeding the energy threshold 
might not produce pairs. To make this determination, we 
again use the Thomson cross section as a reference and 
determine if the photon annihilation takes place by ran¬ 
domly drawing one number S 2 , obtaining by boost¬ 
ing to the center of momentum frame and evaluating if 
> S 2 - If the inequality holds true, the code pro¬ 
ceeds with the pair production calculation. Otherwise, 
it is abandoned and the code returns to step 2. 

Once the photons succeed in producing leptons, the po¬ 
lar scattering angle 0'^ of the newly born e“ is computed 
from the pair annihilation differential cross section as 
given by 

^ 1-^^cos^^'+2(^) b^sin^Oj 

2 \ £1' / (1 — cos^ 

( 6 ) 

(see Jauch & Rohrlich 1980, p.300) where b = 

. A random azimuthal angle (/)(, € [0, 27r) 

is assigned to the e~. Note that Lorentz transformation 
to the zero momentum frame is necessary because equa¬ 
tion m is frame dependent. Utilizing conservation laws, 
the four momenta of the can be determined. 


2.2.3. Process 3: Pair Annihilation/ Photon Production 

The pair annihilation channel is chosen if the random 
particle selection constitutes an e“e+ pair. As with the 
other channels, we first determine the incident angle 9ee 
(subscript ee stands for lepton pair annihilation) between 
the pair by computing the probability distribution of 
scattering as given by 

£^ee(/^e~ ;/^e+; ^ee) ^ElOeefkin (7) 

where fkin as obtained from (Coppi & Blandford 1990) 
is given by: 

fkin = \JPi- + 1^1+ - 131-131+ sin^ 9ee - 2/3e-;5e+ COS Bee- 

( 8 ) 

Here Pe = ^ be. the ratio of lepton speed to the speed of 
light. The code transforms all quantities to the rest frame 
of the electron to calculate the the total cross section aee 
as (Jauch & Rohrlich 1980, p.269): 


( 7 ' + 7 + 4 ) In ( 7 ' + v9^) - P'il' + 3)' 


y(y + 1) 


( 9 ) 


where P' = u'_|_/c, 7 ' = 




i.e. the e"*' speed and 


Lorentz factor respectively in the co-moving frame trav¬ 
eling with the e~. On comparing the tTee/CTT with a ran¬ 
dom number S 3 the code evaluates the occurrence of the 
annihilation event. If the event fails, the code returns to 
step 2 to re-select another pair of particles. Following a 
successful event, the polar scattering angle between ei¬ 
ther of the pair produced photons is generated from the 
differential cross section (from Jauch & Rohrlich 1980, 

p.268) 
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where x = cos0( and d = 7'(1 — P'x). As pointed out 
in the preceding sub-sections, Lorentz transformation to 
the electron frame is necessary as equation (HUD is ex¬ 
pressed in terms of quantities defined in the electron’s co¬ 
moving frame. The random azimuthal angle p'g € [0, 27r) 
is randomly assigned to either photon. Using conserva¬ 
tion laws, the four momenta of the pair produced photons 
can be obtained. 


2.3. Step 3: Back to the lab frame 

At the end of the event, the code transforms the four 
momenta of the particles back to the lab frame by em¬ 
ploying Lorentz transformations. The loop is repeated 
until equilibrium is restored i.e. when the particle num¬ 
bers saturate and distributions become thermal. 

3. RESULTS 

We employ the Monte Carlo code described above to 
study the evolution of the radiation spectrum in a closed 
box containing leptons and photons. The simulations are 
initialized with a Wien radiation spectrum at 10® K and 
a non-equilibrium lepton population, either because lep¬ 
tons and photons are at different temperature or because 
the leptons energy distribution is non-thermal. This is 
expected to mimic a scenario in which the leptons and 
radiation were initially at equilibrium, but the lepton 
population has been brought out of equilibrium by a sud¬ 
den energy release. Such energy release may be due to 
shocks in the fluid (e.g., Rees & Meszaros 1994; Laz- 
zati & Begelman 2010) or by magnetic reconnection in a 
magnetized outflow (e.g. Giannios & Spruit 2006, McK¬ 
inney & Uzdensky 2012). As it will be clear at the end, a 
fundamental parameter that determines the interaction 
between the photons and leptons is the particle ratio, 
i.e., the ratio of photon and lepton number densities. In 
a GRB outflow, such a ratio can be readily estimated. 

Let us call Ek the kinetic energy of the outflow carried 
by particles with non-zero rest mass and E^ the energy 
in electromagnetic radiation. We have: 

E-y NyhVp\^ 

% f 

’^P + ^VlMeV 

By calling rj = E^/{Ey -\- Ek) the radiative efficiency 
of the outflow, and assuming that matter and radiation 
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are coupled in the optically thick region and occupy the 
same volume, equation (El) can be inverted to yield: 

, I (W) 

"'•» 1 (w) 

where the top line is valid for a non-pair enriched fireball 
while the bottom line is for a pair-dominated fireball. All 
values in between are allowed for a partially pair-enriched 
hreball. Note also that we used the convention r 2 = 
r/10^. GRB fireballs are therefore photon-dominated, 
even if highly pair-enriched. 

We here consider two possible values of the parti¬ 
cle ratio. As a representative of pair-enriched plasma, 
we explore the case n.y/niep = 10. A non-enriched 
plasma (or photon-rich plasma) is represented by the ra¬ 
tio n.y/niep = 1000. Note that the latter value is not 
as extreme as the one in equation (11211 . It is, however, 
technically challenging to simulate any higher value of 
the particle ratio. To ensure that the statistics of the 
lepton population is under control, we need to simu¬ 
late at least 1000 irreducible electrons (electrons that 
are not possibly annihilated by a positron). For a parti¬ 
cle ratio rij/niep = 10®, that would require the simula¬ 
tion of 10® photons. We believe that the adopted value 
n-ylnxep = 1000 does capture the characteristics of the 
spectrum emerging from a photon-rich plasma and we 
will discuss the consequences of higher particle ratios in 
Section 4. 

For each particle ratio, we explore different scenarios 
in which the accelerated leptons are either thermal (Laz- 
zati & Begelman 2010) or non-thermal (e.g. Giannios 
2006; Pe’er et al. 2006; Beloborodov 2010) and we con¬ 
sider the possibility of multiple acceleration events, in 
which the leptons are re-energized before the equilibrium 
is reached. Some of these possibilities have been previ¬ 
ously explored, in particular the Comptonization from 
a non-thermal population of electrons (e.g. Pe’er et al. 
2006). We do not consider in this study continuous en¬ 
ergy injection, in which a stationary equilibrium between 
photons and electrons is reached, and for which our code 
is not well suited (e.g. Giannios 2006; Pe’er et al. 2006). 

All simulations are run until equilibrium is attained. 
Here we define equilibrium as the time at which the spec¬ 
tral shape does not change with further collisions and 
the number of photons and leptons saturate. This is 
generally much later than the time at which the total 
energies in leptons and photons approach their asymp¬ 
totic values, since a very small amount of energy can 
make a significant difference in the tails of the distribu¬ 
tion, which are the interesting aspect of the spectrum for 
this study. Our simulations do not have a time stamp, 
since all processes involved are scale free. A time stamp 
can be added upon deciding on a particle and photon 
density, rather than a total number as specified in the 
code. A meaningful comparison with the data can be 
accomplished by considering that a photon in a relativis¬ 
tic outflow with Thomson opacity r scatters-off/collides 
with leptons an average number of times rise — t be¬ 
fore being detected by an observer at infinity (e.g. Pe’er 
et al. 2005). Here we adopt as the Thomson opacity 
of a medium t = f niepXjTds (see Rybicki & Lightman 
1979). It is possible therefore to look at our spectra in 
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Figure 1. Radiation spectrum (upper panel) and leptons’ kinetic 
energy distribution (lower panel) at different simulation stages for 
a photon-rich plasma {N^/N\gp = 1000) with a sudden injection 
of thermal energy in the lepton population (see Sect. [Tl.lll . The 
legend displays the various optical depths at which if energy was 
injected, the corresponding color coded spectrum and distribution 
would be observed. 

the following manner: if a shock or a reconnection event 
dissipates energy in the outflow at a certain optical depth 
r, the spectrum observed at infinity is the one derived 
from our code after Tdiss scatterings per photon. 

3.1. Photon Rich Plasma 

We first explore a photon-rich plasma with n^/ni^p = 
1000. Three simulations are initialized with an out-of- 
equilibrium electron population (there are no positrons 
initially in the plasma) with different initial distributions. 
We inject identical amounts of total kinetic energy K 
in all three cases, raising the average kinetic energy of 
the leptons to 1.365 MeV. This can be considered as a 
mild energy injection and within the equipartition shock- 
acceleration scenario this corresponds to either a mildly- 
relativistic shock or a relativistic shock with a fairly 
low fraction of energy given to electrons (ce <C 1) (e.g. 
Guetta et al 2001). In the hrst simulation, the leptons 
adopt a Maxwell-Jiittner distribution at = 6.5x 10® K. 
In the second case the leptons conform to a Maxwellian 
distribution (at 10® K) which is smoothly connected to 
a non-thermal power law tail Ueij) oc 7 “® ®. The third 
and final simulation explores the scenario of energy dissi¬ 
pation via multiple (10) less energetic injections instead 
of a single intense injection event. 

3.1.1. Thermal Leptons at 6.5 x 10® K 

The lepton population in this case is shocked and then 
thermalizes at 6.5 x 10® K. A similar scenario was ex¬ 
plored analytically and with a simplified Monte Carlo 
code by Lazzati & Begelman (2010). 

The results of the simulation are shown in Figure [1] 
where the evolution of the radiation spectrum and of the 
spectrum of the kinetic energy of the leptons’ population 
are displayed. We first note that the hnal distributions 
(blue curves in both panels) are all thermal, as expected 
for a plasma in equilibrium. Looking at the intermedi¬ 
ate spectra in more detail, we notice that the immediate 
reaction of the radiation spectrum is the formation of a 
high-frequency non-thermal tail, initially appearing as a 
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Figure 2. Evolution of the Band parameters q, 0 and, Ep of 
spectra from the simulation shown in Figure^ The x-axis indicates 
the optical depth of energy injection. 



Figure 3. Fitting of the Band parameters a, 0 and Ep of spectra 
from the simulation shown in Figure fTI at Tdiss = 103. 

new component (for spectra at r^iss = 0 . 001 ) and sub¬ 
sequently forming a continuous tail stemming from the 
thermal photon population {rdiss =5)- At a subsequent 
stage, the low-frequency part of the radiation spectrum is 
also modified, with the spectral peak migrating to higher 
frequencies and causing a flattening of the low-frequency 
component {Tdiss = 75). The hgure shows that the spec¬ 
trum takes a very large number of scatterings for equi¬ 
librium restoration, especially for frequencies lower than 
the peak. For energy dissipation at optical depths up to 
^ 100 a high-frequency non-thermal tail is observed. A 
non thermal low-frequency tail is instead observed even 
for a larger optical depth, up to a few thousand. 

In order to quantify our synthetic transient spectra and 
compare them with observations, we fit them to an an¬ 
alytic model. We adopt the widely used Band function 
(Band et al. 1993) and ht it to the data over a fre¬ 
quency range of three orders of magnitude. Although 
the GRB spectra are in most cases more complex than 
a Band function (e.g. Burgess et al. 2014, Guiriec et 
al. 2011, 2013) this still constitutes a zero-order test 
that any model should pass. We begin by computing 
the mean frequency from our data and select neighbor¬ 
ing frequencies within 1.5 orders of magnitude around 
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Figure 4. Color coded photon spectrum (upper panel) and lep¬ 
tons’ kinetic energy distribution (lower panel) at differe nt stages for 
the photon-rich simulation discussed in Section 13.1.21 The legend 
displays the various optical depths at which if energy was injected, 
the corresponding color coded spectrum and distribution would be 
observed. 

the mean. This data set is binned in frequency and a 
best-fit Band function is obtained by minimizing the . 

Figure [2] shows the evolution of the spectral parame¬ 
ters a, P and Ep for increasing optical depths. We again 
emphasize that this should not be considered as a time 
evolution, since the number of scattering is set by the 
optical depth at which the energy is released in the lep¬ 
tons. A sample fit of the spectrum at Tdiss = 103 to 
the Band function is shown in Figure |3l The hgure rep¬ 
resents a typical case, and shows that the Band model 
hts well the frequencies around the peak but deviations 
are observed for the lowest and highest frequencies. We 
will address this issue further in the discussion. The leg¬ 
end at the top of the hgure shows the Band parameters 
for the ht. An interesting aspect of these simulations is 
that the low-frequency photon index a and the peak fre¬ 
quency are strongly anti-correlated. This is due to the 
fact that it is necessary that the peak frequency shifts to 
higher values for the low-frequency spectrum to change 
from its thermal equilibrium shape. We also note that 
the high-energy slope anticipates the low-energy one, the 
non-thermal features building-up earlier and disappear¬ 
ing faster. We will discuss in more detail these correla¬ 
tions and their implications in Section 4. 

3.1.2. Maxwellian leptons at 10® K with a power law tail 

p = 2.2 

Most models of internal shocks predict the accelera¬ 
tion of non-thermal particles. Gomptonization of seed 
thermal photons by non-thermal leptons has been widely 
studied in different scenarios and under different assump¬ 
tions (e.g. Giannios 2006; Pe’er et al. 2005, 2006). In 
this scenario the shock generates a non-thermal lepton 
distribution characterized by 

N{E)dE oc 7 “^c ?7 (13) 

where 7 is the lepton Lorentz factor and p = 2.2. The 
results of the simulation are shown in Figure 01 where we 
present the evolving radiation spectrum and distribution 
of the kinetic energy of the leptons’ population. We no¬ 
tice that the equilibrium photon and lepton distributions 
(blue curves) are thermal, as expected at equilibrium. 
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Figure 5. Evolution of the Band parameters a, 0 and Ep of 
spectra from the simulation shown in Figure[4] The x-axis indicates 
the optical depth of energy injection. 

We also notice that the spectrum appears non-thermal 
for a wide range of opacities. Initially a prominent high- 
energy power-law tail is developed, for a very small opac¬ 
ity (or Tdiss ~ 0.01). As the injection opacity increases, 
the power-law tail is truncated at progressively lower fre¬ 
quencies, the peak frequency shifts to higher values, and 
a non-thermal tail at low-frequencies develops. The high- 
frequency tail disappears for Tdiss 400, but even larger 
opacities are required to turn the low-frequency tail back 
to the scattering-dominated equilibrium spectrum. We 
fit the Band function to our synthetic spectra and ob¬ 
tain Figure [HI which shows the evolution of the spectral 
parameters a, j3 and Ep for increasing injection optical 
depths. We also notice correlations between the spec¬ 
tral para meters a and the peak frequency, as discussed 
in Section ld.l.ll 

3.1.3. Discrete Multiple Energy Injections 

The presence of multiple minor shocks has been em¬ 
phasized in 2D axisymmetric numerical simulations of 
jets in collapsars (e.g. Lazzati et al. 2009) and seem 
to be an even more common feature in 3D simulations 
(Lopez-Camara et al. 2013). Hence, to provide a more 
realistic scenario for the energy injection we explore lep¬ 
ton heating by multiple energy injections mimicking mul¬ 
tiple shocks instead of a single more powerful one. The 
total energy injected into the lepton population is identi¬ 
cal to the amount injected in the simulations discussed in 
Sections [3.1.H and 13.1.21 However, the energy is divided 
into 10 equal and discrete partitions with each one being 
injected and distributed uniformly among the leptons, 
after every million scatterings. 

The results of the MC simulation are shown in Fig¬ 
ure where the evolving radiation spectrum and the 
spectrum of the kinetic energy of the leptons’ population 
are displayed. In comparison to Figures [1] and 01 two 
differences are apparent for small optical depths. First, 
the high-frequency tail develops much more slowly. Sec¬ 
ondly, the slowly developing tail does not extend to the 
same high energies and in fact, it never approaches the 
MeV mark. Neither of these differences is surprising, 
given that a smaller amount of energy is injected at reg¬ 
ular intervals. The results of the Band function fitting 
are reported in Figure [7] and bring to our attention that 
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Figure 6. Photon spectrum (upper panel) and leptons’ kinetic en¬ 
ergy distribution (lo wer pa nel) at different stages of the simulation 
discussed in Section l.3.1.31 The legend associates the various opti¬ 
cal depths of energy injection with the corresponding color coded 
spectrum and distribution observed. 
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Figure 7. Evolution of the Band parameters a, and Ep of 
spectra from the simulation shown in Figure[6| The x-axis displays 
the opacity at which energy deposition occurred. 



Figure 8. Magnified version of Figure [7| depicting the response of 
the Band function parameters a, /3 and Ep to discrete and multiple 
energy injections, indicated by the broken black vertical lines. The 
x-axis displays the opacity at which energy deposition occurred. 
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like previous other simulations, the high-frequency pho¬ 
ton index /? is the first to respond, and also the first to 
drop just when the a parameter reaches it’s minimum 
value. Another remarkable aspect of the multiple injec¬ 
tion scenario is the immediate reaction of the spectrum 
to new injections, especially for the high-frequency pho¬ 
ton index and the peak frequency (see Figure [5]). 

What is perhaps mostly interesting, rather than the sub¬ 
tle differences among the three scenarios discussed here, 
is the fact the Band parameters of Figures O O and [7] 
show remarkably similar behavior, even though the injec¬ 
tion scenarios are very different. In all three cases, injec¬ 
tion at low optical depth only produces a high-frequency 
power-law tail. Injection at moderate optical depths 
(jdiss 10 — 100) produces a high-frequency power-law 
tail, a shift in the peak frequency, and a non-thermal 
low-frequency tail. Injection at high to very high optical 
depths only results in a non-thermal low-frequency tail 
(see also Section 4 for a discussion). 

3.2. Pair Enriched Plasmas 

In this section we investigate plasmas enriched by 
pairs, by choosing n^/niep = 10. GRB plasmas 
can become pair enriched via energy injection through 
shocks/magnetic dissipation (Rees & Meszaros 2005; 
Meszaros et al. 2002; Pe’er & Waxman 2004) and if the 
peak energy of the resulting distribution exceeds 20 keV 
(Svensson 1982). The generation of pairs is also evident 
from the photon and lepton distributions crossing the 511 
keV mark as shown in the simulations in Sections 13.1.11 
and 13.1.21 We assume, as in the previous scenario, that 
the the pair enriched leptons are impulsively heated by 
injecting equal amounts of kinetic energy AT/IO for the 
first two simulations, albeit with different distribution 
functions (Maxwellian and Maxwellian plus power law). 
The third simulation explores the spectral evolution of a 
pair enriched plasma with an even greater kinetic energy 
injection. The initial photon count of the plasma is 
1.01 X 10®. Being pair-enriched, the total lepton count 
Ne of the plasma is 1.01 x 10^, 

Ne = N^- + 2iVe+e- = 10^ + 10^ (14) 

where N^,- are electrons associated with protons and 
Nf,+e- denotes the number of pairs in the system. 

3.2.1. Maxwellian leptons 

We initiate the simulation with Maxwellian pair- 
enriched leptons that have been impulsively heated to 
10® K, thereby taking the population out of equilibrium 
with the photons. The results of the simulation are dis¬ 
played in Figure [9] with the upper panel depicting the 
photon spectra and the lower panel illustrating the ki¬ 
netic energies of the leptons. Firstly, as observed in the 
section on photon rich plasmas, the final (blue curve) 
spectra is consistent with the equilibrium Wien distri¬ 
bution. For Tdiss ~ 1 a bump is observed to spike near 
the annihilation line along with a power law tail (black 
curve). The lepton distribution also displays a two com¬ 
ponent distribution (black curve in the lower panel). For 
Tdiss 2.3, the power law tail extends farther to high fre¬ 
quencies and merges with the annihilation bump (cyan 
curve). On increasing the injection opacity to around 13, 
the low frequency spectrum flattens, the peak frequency 



Energy (keV) 


Figure 9. Photon spectrum (upper panel) and leptons’ kinetic 
energy distribution (lower panel) at di fferent stages of the pair- 
enriched simulation discussed in Section [3.2.1l The legend displays 
the various optical depths at which if energy was injected, the 
corresponding color coded spectrum and distribution would be ob¬ 
served. 



Figure 10. Evolution of the Band parameters a, /3 and Ep of 
spectra from the simulation shown in Figure|9]for increasing values 
of energy-injection optical depths. 


increases and the annihilation bump merges completely 
with the initial Wien distribution (or the remnant of 
the initial spectrum) creating a non-thermal flattened 
plateau-like feature (yellow curve). The high-frequency 
power law tail returns to the equilibrium Wien spectrum 
much earlier (r^iss < 32) than the non-thermal low fre¬ 
quency tail, which requires about {tcHss 100) to form 
the equilibrium spectrum. We interpret this behavior to 
the inability of the plasma to support a large population 
of pairs. As a consequence the pairs quickly annihilate 
and a large amount of ~ 511 keV photons are injected in 
the plasma. 

The Band parameters obtained by fitting the Band func¬ 
tions to the simulation spectra are plotted in Figure fTOl 
We note that for moderate optical depths, a = —0.75 
and /3 = —1.15 which corresponds to an extremely non- 
thermal spectrum. We also observe from the lower panel 
of Figurefinithat Ep = 20—40 keV. Furthermore, an anti¬ 
correlation is observed between the Band parameters a 
and j3 and between a and Ep. 
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Figure 11. Photon spectrum (upper panel) and leptons’ kinetic 
energy distribution (lower panel) at different stages of the simu¬ 
lation discussed in Section 13.2.21 The legend associates the vari¬ 
ous optical depths of energy injection with the corresponding color 
coded spectrum and distribution observed. 





Figure 12. Evolution of the Band parameters a, /3 and Ep of 
spectra from the simulation shown in Figure m\ The x-axis dis¬ 
plays the opacity at which energy injection occurred. 

3.2.2. Maxwellian leptons at 10® K with a power law tail 

This simulation initializes the lower energy lepton pop¬ 
ulation as thermally distributed at 10® K and a higher 
energy population with a power law tail. However, pair 
enrichment and constraining the injected kinetic energy 
to K/10 lowers the average kinetic energy per lepton in 
comparison to the photon-rich plasmas. As a result the 
leptons are generated according to the distribution 

N{E)dE oei'j-iyPdi^-l) (15) 

where 7 is the lepton’s Lorentz factor and p = 2.2. The 
red curve in lower panel of Figure 1111 displays the ini¬ 
tial kinetic energy distribution of the lepton population. 
Note that the power law tail does not extend to high ener¬ 
gies as the tail in Figure |4] does. The figure also shows the 
evolution of the photon spectra and leptons’ kinetic en¬ 
ergy as equilibrium restoration occurs. For the photons, 
the initial spectra (red curve) and equilibrium spectrum 
(blue curve) fit the Wien distribution. As is expected, 
pair annihilation produces a hump in the vicinity of the 
511 keV region. Meanwhile, the photons forming the ini¬ 



Energy (keV) 


Figure 13. Photon spectrum (upper panel) and leptons’ kinetic 
energy distribution (lower panel) at diff erent stages of the pair- 
enriched simulation discussed in Section 13.2.31 The legend asso¬ 
ciates the various optical depths of energy injection with the cor¬ 
responding color coded particle spectrum and distribution. 

tial Wien spectrum form a power law tail. Similar to the 
previous scenario, at around Tdiss 2 , the power law 
extends to high frequencies and merges with the grow¬ 
ing annihilation hump (cyan curve). We also observe 
a two component distribution in the lepton panel. By 
Tdiss ~ 13, the two component spectrum transforms into 
a broad band flat-plateau like spectrum (yellow curve) 
with the low-frequency spectrum being modified as well. 
The high frequency spectrum of the magenta curve (for 
Tdiss ~ 45) assumes the exponential cut-off of the Wien 
spectrum while the low-frequency tail is still prominent. 
These features make the transient spectra highly non- 
thermal. 

A comparison of Figure[9]with Figure [TT] informs us that 
the spectra of these two scenarios are quite similar. Con¬ 
sequently, a comparison among Figure ITOl and Figure [12] 
also exhibits very similar results - including the anti¬ 
correlations between a and peak frequency and between 
a and /3. 

3.2.3. Maxwellian leptons at 10® K with a power law tail 

p = 2.2 

This section explores the system when pair en¬ 
riched leptons are distributed according the Maxwell- 
Boltzmann distribution at 10® K for lower energies 
whereas the high energy ones form a power-law tail with 
index p = 2 . 2 . 

Similar to the previously discussed cases, the photon 
spectrum fits the Wien spectrum at equilibrium in 
Figure [131 A remarkable difference between Figure [TS] 
and between Figures [9] and [TT] is that the high frequency 
power law tail catches up with the pair-annihilation 
much earlier {rdiss « 0.05) as depicted by the black 
curve. Remnants of the hump are visible in the black 
and cyan curves. Furthermore, for less than 1 scat¬ 
terings, the low frequency tail becomes softer than 
the Wien spectrum (cyan curve). Another important 
non-thermal feature is the broadband nature of the 
flattened spectrum (the yellow curve extends over four 
orders of magnitude in frequency). By about Tdiss ~ 14, 
the truncated high frequency tail approaches the expo¬ 
nential cut-off of the Wien spectrum, whereas the soft 
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^diss 

Figure 14. Evolution of the Band paramet ers a, /3 and Ep of 
spectra from the simulation shown in Figure 1131 with increasing 
energy-injection opacity. 



Figure 15. Evolution of lepton count A^e fo r the simulations in 
Section |3. 2. II (labeled as Thermal) and 13.2.31 Note that th e lep- 
ton count evolution of the simulations discussed in Sections 13.2 II 
and 13.2.21 is indistinguishable. The x-axis displays the opacity at 
which energy injection occurred. 

low frequency tail still persists. 

The best-fit Band function obtained by minimization 
technique, produces highly non-thermal spectral indices 
(a and f3) but the peak frequency as shown in Figure [TT] 
is relatively high for GRBs. The lack of smoothness 
in the a values for moderate optical depths is due to 
the flatness of the photon spectrum as seen from the 
yellow curve in Figure 1131 which occurs in conjunction 
with the transient saturation phase in the lepton count 
(see Figure [la. Figure [15] also displays and compares 
the lepton count for the simulation in Section 13.2.11 (the 
curves labeled as Thermal, which are indistinguishable 
from the pair evolution in Section 13. 2. 2D . Although 
the initial lepton content of the plasmas in the three 
discussed simulations is identical, the plasma with a 
greater kinetic energy injection can sustain pairs for 
larger optical depths leading to a much broader and 
flatter spectrum. For moderate number of scatter¬ 
ings, we obtain a = —1 and /3 = —0.95. Again, an 
anti-correlation is found to exist between the param¬ 
eters a and /3 and also between a and the peak frequency. 



Figure 16. Plot of Band parameters Ep and a for the various sim¬ 
ulations discussed. The solid curves represent photon-rich plasmas 

( ^ = 1000) whereas the broken curves are indicative of pair- 

^'lep 

enriched plasmas where ( = 10). Note the similarity among 

^'lep 

the curves and the exhibited anti-correlation. 



Figure 17. Plot of Band parameters /3 and a for the various sim¬ 
ulations discussed. The solid curves represent photon-rich plasmas 

= 1000) whereas the broken curves are indicative of pair- 

^^lep 

enriched plasmas where ( ^ = 10). Note the complex behavior 

of the curves, especially the evolution of (3. 

4. SUMMARY AND DISCUSSION 

We present Monte Carlo simulations of Compton scat¬ 
tering, pair production, and pair annihilation 
in CRB fireballs subject to mild to moderate internal 
dissipation. We explore cases of photon-rich media - as 
expected in baryonic fireballs - and of pair-dominated 
media. The leptonic component in our simulations is 
initially set out of equilibrium by a sudden injection of 
energy and the spectrum is followed as continuous col¬ 
lisions among photons and leptons restore equilibrium. 

We find that non-thermal spectra arise from transient 
effects. Such spectra could be advected by the expanding 
fireball and released before equilibrium is reached if the 
dissipation takes place at optical depths of up to several 
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hundred. We show that the transient spectra can be rea¬ 
sonably fit by a Band function (Band et al. 1993) within 
a frequency range of 2-3 orders of magnitude around the 
peak and could therefore explain GRB observations. As 
suggested by Lazzati & Begelman (2010), non-thermal 
features can arise even if both the photon and lepton 
distributions are initially thermal, provided that they are 
at different temperatures. As a matter of fact, we find 
that the spectrum emerging from the fireball after a dis¬ 
sipation event at a certain optical depth does not depend 
strongly on the way in which the energy was deposited in 
the leptons. For the photon-rich cases, the first reaction 
of the photon spectrum to a sudden energy injection into 
the leptons is the formation of a high-frequency power- 
law, either because non-thermal leptons are present or 
through the mechanism described in Lazzati & Begel¬ 
man (2010). If the injection happens at moderate opti¬ 
cal depths, the peak frequency of the photon spectrum 
also shifts to higher frequencies and a non-thermal low- 
frequency tail appears. If the energy injection occurs at 
somewhat large optical depths, the high frequency tail 
disappears and the spectrum presents a cutoff just above 
the peak. The low-frequency non-thermal tail is how¬ 
ever very resilient and only if the dissipation takes place 
at very large optical depths, the equilibrium Wien spec¬ 
trum is attained. The pair-enriched simulations show 
a more complex behavior at low optical depths due to 
pair processes, however we still observe the low frequency 
tail’s resilient behavior. We show that this phenomenol¬ 
ogy is rather independent on the details of the energy 
dissipation process and generated lepton distributions: 
non-thermal leptons, high-temperature thermal leptons, 
and multiple discrete injection events all produce simi¬ 
lar spectra. For the case of the pair-enriched simulations 
however, we obtain peak frequencies that are somewhat 
large in the comoving frame (several hundred keV) mak¬ 
ing this scenario less interesting for explaining observed 
burst spectra. However their complex behavior and ex¬ 
treme peak energies offer a tantalizing explanation for the 
rich diversity observed in peak energies of GRBs (Gold¬ 
stein et al 2012) especially when the peak energies > 
MeV. 

The conclusion we can glean from this study is there¬ 
fore that comptonization of advected seed photons by 
sub-photospheric dissipation continues to be a viable 
model to explain the prompt gamma-ray bursts spec¬ 
trum. Agreement is particularly strong when the dis¬ 
sipation occurs at moderate optical depth (of the or¬ 
der of tens) so that both a high- and a low-frequency 
tail are produced. Dissipation at too low optical depth 
would only produce a high-frequency tail, while dissipa¬ 
tion at too large optical depth would only produce a low- 
frequency tail. In a GRB dissipation is likely to occur at 
all optical depths (e.g. Lazzati et al. 2009). The dissipa¬ 
tion events that occur at moderate optical depth would 
therefore be those mostly affecting the spectrum and giv¬ 
ing it its non-thermal appearance. Bursts characterized 
by a Band spectrum over more than three orders of mag¬ 
nitude of frequency remain however challenging for this 
model, and other effects need to be invoked to avoid devi¬ 
ations from the pure power-law behavior at very low and 
high frequencies. Among these effects, some studied in 
the literature are sub-photospheric, radiation mediated 
multiple shocks (Keren & Levinson 2014), line of sight 


effects (Pe’er & Ryde 2011) and high-latitude emissions 
(Deng & Zhang 2014). 

4.1. Spectral correlations 

Besides finding that the overall shape of the partially 
Gomptonized spectra is qualitatively analogous to what 
observed in GRBs, we find that this model predicts the 
existence of two correlations that can be used as a test of 
its validity. We first notice an anti-correlation between 
the low-energy photon index a and the peak frequency. 
The correlation is clearly seen in Figure 1161 where 
results from all simulations are shown simultaneously. 
All simulations start with the same injected photon 
spectrum, the common point in the lower right of the 
diagram. The leptons in all three of the photon-rich 
simulations are energized to identical total kinetic 
energies K albeit different distribution functions. It is 
clear that the evolution of all photon-rich simulations 
is virtually indistinguishable from each other. As 
more and more scatterings occur, the peak frequency 
initially grows and the low-frequency slope flattens. 
At moderate optical depths (~ 100 in all three cases) 
the peak frequency reaches its maximum, the high 
frequency tail disappears (shown in the Figure flTll and 
the low-frequency tail begins to thermalize, dragging the 
peak frequency to slightly lower values. The correlation 
has two branches, a steeper one for r < 100 and a 
flatter one at r > 100. The second branch corresponds, 
however, to spectra without a high-frequency tail and 
is therefore not expected to represent observed GRBs. 
A similar pattern is followed by the pair-enriched cases, 
with the main difference that larger peak frequencies 
are attained along with softer values for a and /3. The 
evolutionary curves for the pair-enriched cases show 
complexity due to the presence of pairs es peciall y at low 
opacities - with the simulation in Section 13.2.31 showing 
a greater amount of variability due to its ability to 
sustain pairs by temporarily balancing the number of 
pair production and annihilation events (see Figure fTSl) . 
In addition to the a — t'pk anti-correlation, we also 
find hints of an anti-correlation between a and /?. 
This correlation is shown in Figure [T3 and is much 
more complex, reflecting the more complex behavior 
of the high-frequency spectrum with respect to the 
low-frequency one. In the case of the high-frequency 
photon-index /3, the way in which the energy is injected 
in the lepton population matters, each simulation 
producing a different track on the graph. 

Comparing these predictions to GRB spectral data is 
not straightforward, since the correlations should not be 
strong in observational data. Adding together data from 
different bursts, the correlations in the observer frame 
would be diluted by the different bulk Lorentz factors of 
bursts and by the diversity of the particle ratio, radia¬ 
tion temperature, and dissipation intensity among busts 
and pulses in a single burst. Still, some degree of cor¬ 
relation has been discussed in the literature, with con¬ 
tradictory conclusions as to its robustness. The a — i^pk 
anti-correlation has been discussed in large burst samples 
(e.g. Amati et al. 2002; Goldstein at al. 2012; Burgess 
et al. 2014). The anti-correlation has been observed 
for some bursts (Zhang et al. 2011), however it is not a 
common feature among GRBs. 
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Photospheric dissipation models have found it difficult 
to reproduce low frequency photon index a ~ — 1 and 
have been unable to explain the GeV emissions (Zhang 
et al 2011). FigurefT^ displays the emission spectra in the 
rest frame of the burst and once Lorentz boosted the pho¬ 
tons forming the high frequency tail reach GeV energies. 
For low/moderate opacities, our simulations have consis¬ 
tently reproduced the low-energy photon index a < 0 as 
shown in Figure [16] thus providing a possible resolution 
for the mentioned issues. Our current model is unable to 
reproduce a < —1.1 for the parameter space explored, 
however additional effects such can modify and further 
soften the low frequency spectra. Analogous studies of 
comptonization effects in GRB outflows have been per¬ 
formed in the past, for example by Giannios (2006) and 
Pe’er et al. (2006). Our work differs from both of these 
previous studies in both content and methodology. Gi¬ 
annios (2006) studied with Monte Carlo techniques the 
formation of the spectrum in magnetized outflows, con¬ 
sidering a particular form of dissipation and assuming 
that the electrons distribution is always thermalized, al¬ 
beit at an evolving temperature. Pe’er et al. (2006), 
instead, used a code that solves the kinetic equations 
for particles and photons, and considere d inje ction of 
non-thermal particles (as in our Section 13.1.21) as well 
as continuous injection of energy in a thermal distribu¬ 
tion. None of these previous studies consider impulsive 
injection of energy in thermal leptons, as discussed here 
or the case of multiple, discrete injection events. In an at¬ 
tempt to keep our results as general as possible we have 
performed the calculations in a static medium, rather 
than in an expanding jet. As long as the opacity at 
which the dissipation occurs is not too large, this should 
not be a major limitation, and the advantage is that 
our results are not limited to a particular prescription 
for the jet radial evolution. In addition, most of the in¬ 
teresting results (the non-thermal spectra) are obtained 
for small and moderate values of the optical depth (or, 
analogously, of the number of scatterings that take place 
before the radiation is released). It should also be noted 
that the assumption of an impulsive acceleration of the 
leptons that does not affect the photon spectra is likely 
not adequate in a highly opaque medium. A final limi¬ 
tation of this study is that only moderate values of the 
particle ratio can be explored. This is an inevitable lim¬ 
itation when both the lepton and photon distributions 
are followed in the scattering process with a Monte Garlo 
technique. If one of the two significantly outnumbers the 
other, a very large number of photons (or leptons) are re¬ 
quired, making the calculation extremely challenging and 
would require parallelizing the code. While performing 
such simulations is important and will eventually become 
possible, we do not anticipate big phenomenological dif¬ 
ferences with respect to what we consider here. Even 
with less electrons, we expect the formation of a high- 
frequency tail (e.g. Lazzati & Begelman 2010), the subse¬ 
quent shift of the peak frequency accompanied by a flat¬ 
tening of the low-frequency photon index, and complete 
thermalization only after many scatterings (i.e., only if 
the dissipation occurs at a very high optical depth). 
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